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Abstract 

Large-scale three-dimensional numerical simulations of the deflagration stage of a thermonu- 
clear supernova explosion show the formation and evolution of a highly convoluted turbulent 
flame in a gravitational field of an expanding carbon-oxygen white dwarf. The fiame dynam- 
ics is dominated by the gravity- induced Rayleigh- Taylor instability that controls the burning 
rate. The thermonuclear deflagration releases enough energy to produce a healthy explosion. 
The turbulent flame, however, leaves large amounts of unburnt and partially burnt material 
near the star center, whereas observations imply these materials only in outer layers. This 
disagreement could be resolved if the deflagration triggers a detonation. 



1 



According to observations and models, many stars that steadily burn their miclear fuel for millions or billions 
of years suddenly end their lives with a powerful explosion that produces a bright object called a supernova. 
A supernova explosion can be powered either by the gravitational energy released during the core collapse 
of a massive; star, or by the nuclear energy released by explosive thermonuclear burning of a star. Here, we 
focus on thermonuclear supernovae that belong to the Type la (SN la) in the observation-based classification 
(1-3). 

Thermonuclear supernovae are produced by explosions of white dwarfs (WD), small dense stars com- 
posed of carbon and oxygen nuclei and detached degenerate electrons {1,3-10). The term degenerate means 
that electrons occupy all possible quantum states below a certain energy. The hydrostatic equilibrium in a 
WD is supported for the most part by the degenerate-electron pressure that does not depend on temperature. 
White dwarfs form at the end of the evolution of stars whose original masses are less than 8 solar masses 
(Mq). a star can lose a large fraction of its material by ejecting outer layers into space at the final stages of 
evolution. The mass of a remaining WD is always less than the Chandrasekhar limit, IAMq, above which 
a hydrostatic equilibrium of degenerate matter is impossible. An isolated carbon-oxygen WD is stable and 
almost inert because its temperature is not high enough to induce any substantial nuclear reactions. This 
isolated dead star can exist almost indefinitely, slowly cooling down as it radiates its energy into space. 
Observations show, however, that more than 50% of all stars are not isolated. They belong to groups of two 
or more stars that orbit a common center of mass. In a close binary system, a WD can increase its own mass 
by accreting material from its companion star. Such systems are considered to be the most probable SN la 
progenitors, even though the exact nature of the companion star and the details of the mass accretion are 
still unclear (1,3-5,7-9). 

When the mass of the WD approaches the Chandrasekhar limit, any small mass increase results in 

a substantial contraction of the star, and the material near its center is compressed. This increases the 
temperature and accelerates thermonuclear reactions near the center. The released energy further increases 
the temperature, thus further accelerating thermonuclear reactions. This process is slowed down by the 
neutrino, convective and conductive cooling. Nevertheless, the temperature in the WD center rises and 
reaches the point where the energy release overwhelms the energy outfiow. In an ordinary non-degenerate 
star, the energy release would be stabilized by a thermal expansion accompanied by the work against gravity. 
In a WD, however, the initial temperature increase does not affect the degenerate-electron pressure, and, 
therefore, does not lead to any substantial expansion that could slow down thermonuclear reactions and 
prevent the runaway process. Eventually, the temperature increases to the level where the thermal and the 
degenerate-electron pressure components become comparable, and the material begins to expand. At that 
time, however, the expansion is already unable to quench the fast thermonuclear burning ignited in the 
center of a WD. The thermonuclear runaway mechanism in degenerate matter was first described in (11) 
and now is a key component of all plausible SN la scenarios (1,4-9). Depending on the scenario, the ignition 
may occur near the center, off center, or in outer layers of the star. Here, we consider the central ignition 
mode. 

Ignition starts a SN la explosion that lasts only a few seconds, but releases ~ lO^'^ ergs, about as much 
energy as the Sun would radiate during 8 billion years. The energy is produced by a network of thermonuclear 
reactions that begins from ^^C and ^^O nuclei and ends in ^^Ni and other iron-group elements. Large amounts 
of intermediate-mass elements, such as Ne, Mg, Si, S, and Ca, are created as well. The main energy-producing 
reactions occur in a thin layer, called a thermonuclear flame, that propagates outwards. At the beginning, 
the flame is laminar and its propagation velocity is controlled by thermal conductivity. As the flame moves 
away from the center, it becomes turbulent and accelerates. Eventually, the burning can undergo a transition 
from a relatively slow, subsonic regime, called deflagration, into a supersonic regime, called detonation, where 
the reaction front is preceded by a shock wave. Most of the energy released during the explosion transforms 
into kinetic and thermal energies of the expanding material. When the sum of kinetic (Ek) and thermal (Et) 
energies exceeds the potential energy of self-gravitation Eg, the star becomes unbound, such that expanding 
material is no longer bound by gravity and will continue to expand indefinitely. 

Thermonuclear reactions that occur during the explosion provide energy for the expansion, but not for 
the luminosity of the expanding gas observed as a SN la. The energy source for the luminosity is the slow 
radioactive decay sequence ^^Ni— >^^Co— >^^Fe. The luminosity reaches its maximum 15 to 20 days after 
the explosion, and then decreases slowly until all the ^^Co decays. The maximum brightness of a SN la is 
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comparable to the brightness of an entire galaxy, and can vary by an order of niagnitnde from one; supernova 
to another. Observations also show that the maximum luminosity of SNe la in visible wavelengths correlates 
with the rate at which the luminosity decreases after the maximum (12-14). This and other approximate 
correlations summarized in (2,8) make it possible to use SNe la as standard candles to measure distances 
and estimate cosmological parameters critical for our understanding of the global evolution of the Universe 
(4,5,7-10,15-19). 

The importance of SNe la as standard candles grows as observational techniques improve and estimations 
of cosmological parameters become more accurate. We are approaching the limit, however, where the reliance 
on empirical correlations becomes the main source of uncertainty. Even though the correlation between the 
maximum luminosity and rate of decrease of the luminosity for SNe la has a theoretical explanation (20) 
based on a one-dimensional model, this is at best only a first approximation that does not take into account 
a detailed explosion mechanism and its possible variations. The only way to solve this problem is to study 
details of supernova explosions using multidimensional numerical simulations. 

One-dimensional (ID) numerical models have been extensively used to test general ideas about possible 
explosion mechanisms (21-24,4,6). Delayed-detonation models (25-32), that postulate a deflagration-to- 
detonation transition (DDT) at some stage of the thermonuclear explosion, are most successful in reproducing 
observed characteristics of SNe la. Many important details, however, including the mechanism of DDT, are 
still unknown because SN la explosions are intrinsically three-dimensional (3D) phenomena. 

Only a full-scale 3D numerical model may be expected to reproduce all key features of the explosion that 
involves propagation of a turbulent thermonuclear flame in the gravitational field of a white dwarf. Building 
such a model is a complicated interdisciplinary problem on the leading edge of astrophysics, nuclear physics, 
combustion physics, and computational physics. Full-scale 3D numerical simulations of thermonuclear super- 
nova explosions have become a reality during the last few years (33-35), in great part owing to the progress 
in computational technology. Here, we describe a 3D numerical model, present numerical results for the 
deflagration stage of the explosion, and discuss implications of the results for astrophysics of SNe la. 

Physical and Numerical Model. The numerical model described in more detail in (36) is based on 

reactive Euler equations that represent mass, momentum, and energy conservation laws for an inviscid fluid. 
The thermodynamic properties of the fluid are defined by the equation of state of degenerate matter, which 
is well known from basic theory and includes contributions from ideal Fermi-Dirac electrons and positrons, 
equilibrium Planck radiation, and ideal ions. Reactions involved in the thermonuclear burning of the carbon- 
oxygen mixture are described by a simplified four-equation kinetic scheme (25,33). This kinetic scheme is 
coupled with a flame-capturing technique (37,33) that introduces an additional partial differential equation 
and ensures thermonuclear flame propagation at a prescribed speed S. The resulting set of equations is 
integrated on a Cartesian adaptive mesh using an explicit, second-order, Godunov-type numerical scheme 
(37,38). The mesh is dynamically refined around shock waves, flame fronts, and in regions of steep gradients 
of density, pressure, composition, and tangential velocity. The computational cell size dx varies within 
predefined limits dxmin and dx^ax- 

The large-scale simulations described here do not resolve the physical thickness of a flame that differs 
from the white dwarf radius RwD by up to 12 orders of magnitude. Therefore, the flame speed must be 
provided by an additional subgrid model that takes into account physical processes at scales smaller than 
the computational cell size. We deflne the flame speed S as 

S = ma^(Si,St), (1) 

where the steady-state laminar flame speed Si in carbon-oxygen degenerate matter is a known function of 
temperature, density, and composition (39,40). The speed St of a quasi-steady-state turbulent flame driven 
by the gravity- induced Rayleigh- Taylor (RT) instability depends on the gravitational acceleration g and the 
length scale L (37,41): 

St ^ (2) 

where A = (po — pi)/(po + Pi) is the Atwood number, and po and pi are the densities ahead and behind the 
flame front, respectively. The driving scale L is set equal to one or two computational cell sizes dx. 

Equation (2) is based on the two main properties of a turbulent flame (37,41): self-similarity of the flame 
structure and self-regulation of the flame speed. Self-similarity means that the 3D distortions of the flame 
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surface at difFcrcnt scales arc similar. Sclf-rcgulation means that changing the flame speed at small scales does 
not affect the flame speed at larger scales. This occurs because a higher flame speed at small scales causes 
small flame wrinkles to burn out, thus decreasing the flame surface. The resulting burning rate, defined as a 
product of the flame speed at small scales and the flame surface, does not change. This subgrid model makes 
it possible to reproduce the correct flame propagation in numerical simulations while explicitly resolving 
only the large-scale flame structure. If the resolved flame structure is self-similar and self-regulating, and 
behaves according to the Eq. (2), the subgrid model just extends this behavior to unresolved small scales. 
Shifting the boundary between resolved and unresolved scales by changing the numerical resolution should 
not affect the turbulent flame propagation. The solution obtained should then be independent of numerical 
resohition and on the exact value of St. Numerical convergence tests described in (36) show that for the high 
{dxmin = 2.6 X 10^ Cm) and medium [dxmin = 5.2 x 10^ cm) numerical resolutions the 3D structure of the 
turbulent flame was resolved in sufficient detail to reproduce key dynamic properties of the flame included 
in the subgrid model. Changing L in Eq. (2) from dx to 2dx has only a minor effect on the converged 
solution. The results are self-consistent, reasonably accurate, and can be used to analyze explosion scenarios 
for thermonuclear supernovae. 

Thermonuclear deflagration in carbon-oxygen white dwarf of Chandrasekhar mass. The initial 
conditions for the simulations (see (36) for more details) were set up for a Chandrasekhar-mass WD in 
hydrostatic equilibrium with the initial radius Rwd = 2 x lO^crn, the initial central density = 2 x 

10'^ g/cw? , the uniform initial temperature T = 10^ K, and the uniform initial composition with equal mass 
fractions of ^^C and ^^O nuclei. The 3D computational mesh extended from the WD center x = y = z = 
to x = y = z = 2.6RwD- Thus, we model one octant of the WD assuming mirror symmetry along the a; = 0, 
y = and z = Q planes. The burning was initiated at the center of WD by filling a small spherical region at 
r < 0.015RwD with hot reaction products without disturbing the hydrostatic equilibrium. 

The development of the explosion for the high-resolution case is shown in Fig. 1 by a series of 3D 
snapshots of the fiame and WD surfaces (see also Movies SI and S2). At the beginning, the initially spherical 
fiame ignited in the center of WD propagates outwards with the normal laminar flame speed Si . As it moves 
away from the center, the gravitational acceleration g and the Atwood number A increase. This increases the 
amplitude and rate of development of the RT instability controlled by g and A. Due to the RT instability, 
small perturbations of the flame surface grow and form a few plumes that have characteristic mushroom 
shapes. The turbulent flame speed St also increases with g and A according to Eq.(2), and eventually 
dominates Si in Eq.(l). The flame plumes continue to grow, due partially to the flame propagation and 
partially to gravitational forces that cause the hot, burnt, low-density material inside the plumes to rise 
towards the WD surface. The same gravitational forces also pull the cold, high-density unburnt material 
between the plumes down towards the center. The resulting shear flows along the flame surface are unstable 
(Kelvin-Helmholtz (KH) instability) and quickly develop vortices. These vortices further distort the flame 
surface, and also contribute their energy into the turbulent cascade that creates turbulent motions at smaller 
scales, down to a few dx. 

When the original fiame plumes grow large enough, secondary RT instabilities develop on their surface, 
thus producing the next level of "mushrooms" that also grow and may become subject to the RT instability 
at a smaller scale, etc. These smaller gravity-induced mushrooms interact with the turbulence created by 
the previous generation of larger flame plumes, and also produce some turbulence themselves through the 
KH instability. The resulting complicated turbulent flame surface is shown in Fig. 1. 

As the turbulent flame develops, the energy released by the thermonuclear burning causes the WD to 
expand. The expansion accelerates and becomes nonuniform as the rising plumes approach the star surface. 
We continued the simulations until the star surface reached the computational domain boundary. By that 
time, the radius of the expanding star increased by about a factor of 2.6, the outer layers accelerated to about 
1.2 X 10^ cm/s, and the density of unburnt material near the star center decreased to about 5 x 10^ g/cm^. 
The area around the center still contains a significant amount of unburnt material that sinks at 10^ cm/s 
towards the center between large flame phimcis. The velocity of the large flame plumes is essentially zero 
relative to the expanding matter, that is, the plumes have practically stopped rising. This effect of freezing 
of the RT-instability on large scales due to expansion is also related to freezing of large-scale turbulence (37), 
and contributes to the burning-rate decrease after 1.5 s. 
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Effects of initial composition and background turbulence. Before we eompare the computed super- 
nova explosion with astronomical observations, we need to know the extent to which the results are affected 
by uncertainties in initial model parameters. We examined two properties of the WD that are not well 
defined and arc expected to have the most pronounced effect on the explosion, the initial composition and 
the level of background turbulence. 

In the base case described above, we assumed a uniform composition with the carbon mass fraction 
Xc = 0.5. A WD can also have a core partially depleted of carbon (42). To take this into accoimt, we 
calculated the explosion for Xc = 0.25 in the core and Xc = 0.5 in the outer layers. The core with the 
radius 0.25Rwd contained about 25% of the total mass of the WD. The lower C concentration resulted in 
a slightly slower explosion as shown in Fig. 2 by the total energy Etot, the kinetic energy Ek, the released 
nuclear energy En, and the burnt mass fraction as functions of time. Compared to the base case, it 
takes about 0.2 seconds longer for the slower explosion to release the same amount of energy. There is no 
significant difference in the turbulent flame structure. 

The background turbulence is created by intense convective flows that appear when the accreting WD 
approaches the Chandrasekhar limit and quickly contracts just before ignition. A flame propagating through 
intense turbulent flow can be accelerated to the velocity of turbulent motions, Vt- Two-dimensional simu- 
lations (42) show that Vt estimated as a differential velocity between adjoining eddies with sizes ~ 10^ cm 
can reach 2.4 x 10^ cm/s in central parts of a WD. To model this effect, we replaced Eq.(l) by 

S^max{Si,St,Vt) (3) 

and assumed that Vt has a constant value V^ for r < 0.33Riyd, and linearly decreases to as r approaches 
RwD- The simulations were performed for = 3 x 10^ cm/s. The results (Fig. 2) show that the back- 
ground turbulence significantly increases the energy-release rate and the expansion at the initial stage of the 
explosion. During the first 0.45 s, the explosion with the background turbulence releases the same amount of 
energy as the base case during 0.8 s. The material burns faster because the laminar flame speed (~ 10^ cm/s), 
which played an important role at the initial stage for the base case, was replaced by a higher value of Vt 
according to Eq.(3). The initial flame speed, however, also has a long-lasting effect because the initial stage 
of explosion is critical for the development of the turbulent flame surface. A flame that initially propagates 
with a higher speed competes with the RT instability more efficiently and develops fewer wrinkles on resolved 
scales. This results in a smaller surface area of the fiame. As the flame propagates outwards, Si increases for 
both cases. It becomes higher than Si for the base case, and higher than Vt for the case with the background 
turbulence. St is approximately the same for both cases, but the explosion with the background turbulence 
begins to slow down when the burnt mass fraction reaches 0.15 because of a smaller flame surface area, 
and eventually releases much less energy than the base case. 

If we increase to the unrealistically high value 10^ cm/s, the background turbulence dominates even 
far from the center and the RT instability does not develop at all. The flame surface is almost spherical until 
the end of the simulation, and there is no unburnt material below the flame surface. The energy release and 
expansion are fast (Fig. 2) until the flame reaches low-density layers where the C burning almost stops. The 
flnal values of Etot, Ek, En, and /(, are higher than for V^ = 3 x 10^ cm/s, but lower than for the base case 
in which the highly convoluted turbulent flame continues to burn high-density material in central parts of 
the WD. 

The uncertainty in the maximum energy released by the explosion can be estimated as the difference 
between the base case and the case with V^ = 3x 10^ cm/s (Fig. 2). Even though the burning continues for 
the base case, it begins to slow down and will stop when the density of the unburnt material near the center 
drops below ~ 10^ g/cm^. We expect that by that time, the burnt mass fraction will increase to 0.6-0.7, 
which is 1.7-2 times higher than for the case with the background turbulence. The final released energies 
for these two cases wifl also differ by a factor of 1.7-2 {En 0.87 x 10^^ ergs for 1/^° = 3 x lO'^ cm/s, and 
En ^ (1.5 - 1.7) X 10^^ ergs for the base case). 

Comparison of simulations and observations. Now, after ensuring numerical convergence and ex- 
amining the sensitivity to initial conditions, we can try to compare results of simulations to astronomical 
observations. One important parameter that can be directly compared to observations is the kinetic en- 
ergy of the expanding material. This energy can be calculated from expansion velocities measured using 
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the Dopplor shift of spectral lines. The typical kinetic energy for SN la is (1 — 1.5) x 10"''^ ergs (6). In 
simulations, this energy corresponds to E/. at infinity, when all the nuclear energy released transforms into 
kinetic energy of the expanding material and the work against gravity: Ek = En + £^°oi- The negative initial 
value E^^^ — —0.5 x 10'"'^ represents the binding energy of the star. For our base case, the estimated final 
value Ek = En — 0.5 x 10^^ ~ (1.0 — 1.2) x 10^^ ergs is at the lower end of the typical range. The case 
with the realistic level of background turbulence, = 3 x 10^ cm/s, produces a weak explosion, with final 
Ek — 0.37 X 10^^ ergs. The explosion is also not strong enough {E^ = 0.6 x 10^"'^ ergs) for the high level of 
background turbulence, = 10^ cm/s. For all the cases shown in Fig. 2, the final value of Etot is positive, 
that is, the explosion releases enough energy to unbind the star. The base-case 3D explosion is much more 
energetic than explosions resulted from two-dimensional simulations (43,26,34) because the flame does not 
develop enough surface area in two dimensions as discussed in (44). 

The key feature of the simulations is the highly convoluted turbulent flame surface that allows extensive 
intcrpcnetration of burnt and unburnt materials. The angle- averaged burnt mass fraction at 1.9 s is less that 
1 for any distance r from the WD center (Fig. 3). In particular, 80 — 90% of the material near the WD center 
is unburnt. This material continues to burn, but it will not burn out completely as long as convective flows 
supply fresh unburnt material from outer layers. As the WD continues to expand, the density decreases. The 
deflagration begins to produce intermediate-mass elements when the density of unburnt material becomes 
lower than ~ 5 x lO'' g/cm^, and the burning stops when the density drops below ~ 10^ g/cm^. This 
means that the final cjccta produced by the 3D deflagration model will contain the unburnt material and 
intermediate-mass elements at any distance from the center. 

The unburnt carbon and oxygen that remain between the flame plumes and intermediate- mass elements 
that form at low densities at different r should produce spectral signatures in a wide range of expansion 
velocities, including low velocities close to zero. Analyses of SN la spectra, however, imply C (45-47) and 
O (48) only at high velocities, as would be produced by the acceleration of expanding outer layers. For 
intermediate-mass elements, minimum observed velocities are lower (2), but large enough (^ 10,000 km/s 
for Si) to rule out the presence of these elements near the WD center. In the simulations, the low-velocity 
unburnt material was not present only for the case with an unrealistically high level of background turbulence 
that produced a spherical flame and resulted in a weak explosion. 

The dynamics of the 3D deflagration described here is in agreement with the preliminary simulations 
(33) performed with the same model. These results were independently conflrmed by Reinecke et al. (34,35) 
who used different subgrid model, nuclear kinetics scheme and initial conditions, and obtained similar highly 
convoluted flames with unburnt material near the WD center. 

We thus conclude that the deflagration model of a SN la explosion is incomplete. The most natural 
solution to this problem that would make the results consistent with observations would be to assume that 
the turbulent flame triggers a detonation. A thermonuclear detonation wave could propagate through the 
WD with velocities ~ 10^ cm/s (49,50) and would quickly burn all the material near the center leaving only 
the low-density outer layers unburnt. For the density below 5 x 10^ g/cm'', a detonation would produce 
intermediate-mass elements (25) that are observed in spectra of SNe la. A detonation would also partially 
smooth out composition inhomogeneities that are predicted by the deflagration model and may be incompat- 
ible with observations (51). Remaining asymmetries may account for a weak polarization recently detected 
in SN la spectra (52,53). 

One-dimensional (25,28-32) and two-dimensional (26,27) delayed-detonation models were the most suc- 
cessful in explaining observable characteristics of SNe la. These models, however, use the time for detonation 
initiation as a free parameter because the DDT problem is intrinsically 3D and still unsolved. A large-scale 
3D model also cannot reproduce DDT phenomena that involve physical processes occurring on small unre- 
solved scales. One approach to solving this problem is to study, in much more detail, the types of reacting 
flows created by 3D deflagrations and look for situations that create the right types of "hot spots" that we 
know (54) are the sources of detonation initiation. 
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FIGURES 

Fig. 1. Development of thermonuclear deflagration in carbon-oxygen white dwarf. The gray surface shows 
the turbulent thermonuclear flame. The color scale shows the radial velocity of unburnt material scaled by 
1000 km/s. Distances are scaled by the computational domain size Xmax = 5.35 x 10* cm. Numbers in the 
frame corners are time in seconds (left top), burnt mass fraction fb (right top), and flame surface area scaled 
by a;^„3. (right bottom). High resolution {dxmin = 2.6 x 10^ cm). L = 2dx. 
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Fig. 2. Total {Etot = Ek + Et~ Eg), kinetic {Ek), released nuclear (En) energies and the burnt mass fraction 
fb as functions of time for the base case without background turbulence (1), with background turbulence 
V^P = 3 X 10^ cm/s (2), and = 10^ cm/s (3). Dotted lines correspond to reduced carbon concentration 
Xc = 0.25 in the core. Energy units are 10^° ergs. Medium resolution {dxmin = 5.2 x 10^ cm). L = dx. 




Fig. 3. Aiiglc-avcraged local mass fraction of burnt material (dotted line) and total burnt mass fraction 
inside a spherical region of radius r (solid line) as functions of distance r from the WD center. The distance 
is scaled by the computational domain size Xmax = 5.35 x 10® cm. Time 1.9 s corresponds to the last frame 
of Fig. 1. 
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Supporting Online Material 

V. N. Gamezo, A. M. Khokhlov, E. S. Oran, A. Y. Chtchelkanova, R. O. Rosenberg. Explosions of Ther- 
monuclear Supernovae: Physical Implications of FuU-Scale, Three-Dimensional Simulations. 

MATERIALS AND METHODS 

NRL Supernova Model 

Fluid Dynamics and Equation of State 

The numerical model is based on reactive Euler equations 

V-(pU) , 

V • (pUU) -VP + pg, 
V-{V{E + P))+pV-g + pq , 

where p, E = Ei+pU^ /2, Ei, U,g and q are the mass density, energy density, internal energy density, velocity 
of matter, gravitational acceleration, and nuclear energy release rate per unit mass, respectively. These 
equations describe mass, momentum, and energy conservation laws for an inviscid fluid. The thermodynamic 
properties of the fluid are defined by the equation of state of degenerate matter, which is well described by 
basic theory and includes contributions from ideal Fermi-Dirac electrons and positrons, equilibrium Planck 
radiation, and ideal ions. Pressure P = P{p, Ei, Yg, Yi) and temperature T = T{p, Ei, Y^, Y,) are determined 
by the equation of state as functions of p, Ei, the electron mole fraction Y^, and the mean mole fraction 
of ions Yi. The relation between P and p is close to that in a polytropic gas with 7 varying from 4/3 to 
5/3. This equation of state is valid for the thermodynamic parameters and compositions expected in the 
computations, ranging from relatively cold highly degenerate carbon-oxygen matter to partially degenerate 
hot products of thermonuclear reactions. 

Nuclear Kinetics and Flame Propagation 

Thermonuclear reactions involved in the thermonuclear burning of the carbon-oxygen mixture (S1-S3) 
can be separated into three consecutive stages responsible for the energy release. First, the ^^C + ^^C reaction 
leads to the consumption of C and formation of mostly Ne, Mg, protons, and a-particles. Then begins the 
nuclear statistical quasi-equilibrium (NSQE) stage, during which O burns out and Si-group (intermediate 
mass) elements arc formed. Finally, Si-group elements are converted into the Fe-group elements and the 
nuclear statistical equilibrium (NSE) sets in. The reaction time scales associated with these stages strongly 
depend on temperature and density and may differ from one another by several orders of magnitude {S4-S8). 

The full nuclear reaction network includes hundreds of species that participate in thousands of reactions. 
Integration of this full network is too time-consuming to be used in multidimensional numerical models. 
Therefore, we used a simplified four-equation kinetic scheme {S9,S10) that describes all major stages of 
carbon burning. This kinetic scheme is coupled with a flame-capturing algorithm {S11,S10) that ensures 
thermonuclear flame propagation with a prescribed speed S. Flame capturing is needed because we cannot 
resolve in large-scale simulations the physical thickness of a laminar flame that differs from the white dwarf 
radius Rwd by up to 12 orders of magnitude. 

Numerical Method 

The fluid dynamic equations, coupled to the nuclear reaction mechanism and the flame-capturing 
algorithm are integrated using an explicit, second-order, Godunov-type, adaptive-mesh-refinement code 
{S11,S12). A Riemann solver is used to evaluate fluxes at cell interfaces. The computational mesh is 



dp 

dt 

dpV 

dE 
'dt 
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comprised of cubic cells of various sizes that are organized in a fully threaded tree (FTT) {S12). The 
FTT-based parallel adaptive mesh refinement algorithm dynamically adjusts cell sizes in accordance with 
changing physical conditions in the vicinity of each cell. Here, the mesh was refined around shock waves, 
flame fronts, and in regions of steep gradients of density, pressure, composition, and tangential velocity. The 
cell size dx varies within predefined limits dxmin and dxmax in such a way that neighboring cell sizes can be 
the same or diflFer by a factor of two. The code has been extensively tested and used in various combustion 
problems involving shocks, flames, turbulence, and their interactions {S13,S14, and references therein) and 
in astrophysics (Si 5). 

Initial Conditions 

The initial conditions for the simulations were set up for a Chandrasekhar-mass white dwarf (WD) 
in hydrostatic equilibrium with the initial radius Rwd = 2 x lO^cm, the initial central density Pc = 2 x 
10^ g/cm^, the uniform initial temperature T = 10^ K, and the uniform initial composition with equal 
mass fractions of ^^C and ^^O nuclei. Starting from the central pressure P{pc)- the equations of hydrostatic 
equilibrium, dP/dr = —GMp/r^ and dM/dr = Awpr^, were integrated outward until P = was reached 
(here G is the gravitational constant, and M is the mass of the material inside a sphere of radius r). The 
resulting WD configuration was interpolated onto a 3D mesh extended from the WD center x = y = z = 
to X = y = z = 2.QRiYD- Thus, we model one octant of the WD assuming mirror symmetry along the a; = 0, 
y = and 2; = planes. The burning was initiated at the center of WD by filling a small spherical region at 
r < O.OlbRwD with hot reaction products without disturbing the hydrostatic equilibrium. 

Because a Chandrasekhar-mass WD is close to a collapse threshold, its gravitational equilibrium is 
sensitive to the discretization errors that appear when the spherical body is mapped into a Cartesian mesh. 
To minimize these errors, the mesh is initially refined to the finest level near the WD center (r < OAR'iyd)- 
During the simulations, we keep the mesh unchanged until the flame reaches the flne grid boundary. Then 
the adaptive mesh refinement algorithm is turned on. 

Subgrid Model for Flame Speed 

The speed of a laminar thermonuclear flame in a WD is defined by reaction rates and transport prop- 
erties of the material and governed by the same laws that describe laminar flame structure in terrestrial 
chemical systems {S16-S18). The only substantial difl[erence is that transport properties of degenerate mat- 
ter are dominated by the electron heat conduction at high densities, and by both electron and photon heat 
conduction at low densities. The steady-state laminar flame speed Si in carbon-oxygen degenerate matter 
is a known function of temperature, density, and composition {S19,S20), but the flame can be laminar only 
near the center of a WD. Away from the center, the flame is turbulent and propagates with a higher effective 
speed {S21,S11,S22,S23). 

In our simulations, the turbulent flame speed St is deflned by a subgrid model that takes into account 

physical processes at scales smaller than the computational cell size. The subgrid model assumes that burning 
on small unresolved scales is driven by the gravity-induced Rayleigh- Taylor (RT) instability. This instability 
distorts the flame surface at multiple scales and generates turbulent motions in the surrounding fluid. The 
turbulent energy propagates from large to small scales, thus further disturbing the flame surface on small 
scales. A developed turbulent flame is statistically steady-state and forms a dynamic hierarchical self-similar 
3D structure where the flame surface is distorted at multiple scales. The flame distortions increase the flame 
surface and, therefore, the burning rate. The larger the scale, the higher the burning rate at this scale due 
to the increase of the flame surface resulting from distortions at smaller scales. The burning rate defines 
the turbulent flame speed St for any given length scale. This turbulent flame structure was analyzed in 
3D numerical simulations {S11,S24) performed for the thermonuclear burning of carbon-oxygen degenerate 
matter in a uniform gravitational field. It was found that a turbulent fiame in a vertical column of width L 
becomes quasi-steady-state and propagates with the speed 

St =i 0.5^/A^ {SI) 

independent of the laminar speed Si, where A = {po — pi)/{po + Pi) is the Atwood number, and po and 
pi are the densities ahead and behind the flame front, respectively. We used this result in the subgrid 
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model assuming that at L « Rwd burning can be considered as locally steady state (S25,S26,S11). The 
driving scale L was set equal to the computational cell size dx. The resulting flame velocity S used by the 
flame-capturing algorithm was calculated as 

S = m&x{Si,St) {S2) 

Numerical Convergence 

Numerical convergence tests, which are critical for establishing the validity of any numerical simulations, 
are performed to ensure that the solution obtained is not significantly affected by the computational cell 
size dx. In our case, numerical convergence is closely related to the physics of the subgrid model that 
defines the turbulent flame speed in Eq.(Sl) and assumes that burning can be considered as locally steady 
state. The steady-state assumption is not valid for length scales affected by the spherical geometry and 
expansion. These scales need to be explicitly resolved, but this will not necessarily ensure that the solution 
is independent of dx. The resolved flame properties should also correspond to the flame properties built into 
the subgrid model. 

The subgrid model is based on the two main properties of a turbulent flame {S11,S24): self-similarity 
of the flame structure and self-regulation of the flame speed. Self-similarity means that the 3D distortions 
of the flame surface at different scales are similar. Self-regulation means that changing the flame speed 
at small scales does not affect the flame speed at larger scales. This occurs because a higher flame speed 
at small scales causes small flame wrinkles to burn out, thus decreasing the flame surface. The resulting 
burning rate, deflned as a product of the flame speed at small scales and the flame surface, does not change. 
This subgrid model makes it possible to reproduce the correct flame propagation in numerical simulations 
while explicitly resolving only the large-scale flame structure. If the resolved flame structure is self-similar 
and self-regulating, and behaves according to the Eq.(Sl), the subgrid model just extends this behavior 
to unresolved small scales. Shifting the boundary between resolved and unresolved scales by changing the 
numerical resolution should not affect the turbulent flame propagation. The solution obtained should then 
be independent of numerical resolution and on the exact value of St. 

We performed a series of computations with three resolutions: low [dXmin = 10-5 x 10^ cm), medium 
{dx„iin = 5.2 X 10''' cm), and high {dx„iin = 2.6 x 10^ cm). The maximum computational cell size dx„iax = 
41.8 X 10^ cm and the computational domain size Xmax = 5.35 x 10* cm were kept constant. The development 
of the explosion for the high-resolution case is shown in Fig. 1 and described in the main text. The simulation 
with medium resolution showed a similar pattern of development, though it produced fewer small-scale 
wrinkles on the flame surface. The low-resolution simulation showed only the largest flame plumes. The 
surface of these plumes remained smooth and did not develop any significant instabilities. These three cases 
are compared in Fig. SI, which shows the total energy Eioi, the kinetic energy Ek, the released nuclear 
energy and the burnt mass fraction /b as functions of time for the developing explosions. The curves 
for high and medium resolutions practically coincide for £^tot) En, and /(,. The kinetic energy Ek shows the 
largest difference between the high-resolution and the medium-resolution cases, and this does not exceed 
10%. Therefore, the solution is practically converged. 

In order to test the self-regulating properties of the flame, we increased L in Eq.(Sl) by a factor of 2 
(L = 2dx), and repeated the three simulations with the same high, medium, and low numerical resolutions. 
The key energies and for these simulations are compared with previous three cases in Fig. S2. For 
the low-resolution case, the self-regulation mechanism on resolved scales does not work, St and key energies 
increase by the same factor \/2. For the medium and high resolutions, the self-regulation appears on resolved 
scales and reduces the maximum difference between the energies calculated for different St to 15% and 6%, 
respectively. These tests show that the numerical convergence was achieved in the high- and medium- 
resolution simulations. The 3D structure of the turbulent flame was resolved in sufficient detail to reproduce 
key dynamic properties of the flame included in the subgrid model. The results are self-consistent and 
reasonably accurate, and can be used to analyze explosion scenarios for thermonuclear supernovae. 
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SUPPORTING FIGURES 




Fig. SI. Total (Etot), kinetic (Ek), released nuclear (En) energies and the burnt mass fraction fb as 
functions of time. Solid, dashed, and dotted lines correspond to high, medium, and low numerical resolution, 
respectively. Ek is calculated from local fluid velocities and densities as a total for all computational cells. 
The conservation law Ek + Et — En — Eg = const is satisfied with accuracy better than 1%. 
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Fig. S2. Total (Etot), kinetic (Ek), released 
nuclear energies and the burnt mass 

fraction as functions of time for high (a), 
medium (b) and low (c) numerical resolutions, 
L = dx (solid lines) and L = 2dx (dotted 
hues). 
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SUPPORTING TABLES 



Table SI. Parameters displayed in Movie 2. 



Parameters Definition 

step, fburn timestep, mass fraction of burnt material 

time, rhomax time (s), maximum mass density in computational domain (g/cm^) 

fsurf, fsurfl fiame surface area (^), fsurfl = fsurf/{A'Krfmax'^) 

rfmin, rfmax min, max distance from star center to flame surface (^) 

ssurf, ssurfl star surface area (-^), ssurfl = ssurf /{An rsmax"^) 

rsmin, rsmax min, max distance from star center to star surface (^) 

vrsmin, vrsmax min, max radial flow velocity at star surface (cm/s) 

vrmin, vrmax min, max radial flow velocity in computational domain (cm/s) 

etherm thermal energy (^) 

ekin kinetic energy {^) 

egrav gravitational enegry (^) 

etot total energy = etherm + ekin - egrav (^) 

enuc nuclear energy (^) 



(1) scaled by xl^^^. 

(^) scaled by Xmax = 5.35 x 10* cm. 

(^) scaled by 10^°ergs 
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